Mon. Not. R. Astron. Soc. 000,[Tf?? (2012) Printed 5 December 2012 (MN KTeX style file v2.2) 



3D-PDR: A new three-dimensional astrochemistry code for treating 
Photodissociation Regions 
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ABSTRACT 

Photodissociation regions (PDRs) define the transition zone between an ionized and a dark 
molecular region. They consist of neutral gas which interacts with far-ultraviolet radiation 
and are characterized by strong infrared line emission. Various numerical codes treating one- 
dimensional PDRs have been developed in the past, simulating the complexity of chemical 
reactions occurring and providing a better understanding of the structure of a PDR. In this 
paper we present the three-dimensional code, 3D-PDR, which can treat PDRs of arbitrary 
density distribution. The code solves the chemistry and the thermal balance self-consistently 
within a given three-dimensional cloud. It calculates the total heating and cooling functions at 
any point in a given PDR by adopting an escape probability method. It uses a HEALPix-based 
ray-tracing scheme to evaluate the attenuation of the far-ultraviolet radiation in the PDR and 
the propagation of the far-infrared/submm line emission out of the PDR. We present bench- 
marking results and apply 3D-PDR to i) a uniform-density spherical cloud interacting with a 
plane-parallel external radiation field, ii) a uniform-density spherical cloud interacting with a 
two-component external radiation field, and iii) a cometary globule interacting with a plane- 
parallel external radiation field. We find that the code is able to reproduce the benchmarking 
results of various other one-dimensional numerical codes treating PDRs. We also find that 
the accurate treatment of the radiation field in the fully three-dimensional treatment of PDRs 
can in some cases leads to different results when compared to a standard one-dimensional 
treatment. 

Key words: Keywords - methods: numerical; radiative transfer; (ISM:) H ii regions; ISM: 
abundances; astrochemistry. 



1 INTRODUCTION 

Photodissociation regions (PDRs; also known as Photon Domi- 
nated Regions) are ubiquitously present in the interstellar medium 
(ISM), and consist of predominantly neutral gas and dust illumi- 
nated by far-ultraviolet (FUV) radiation (6 < hv < 13.6 eV). 
Studies of PDRs allow us to understand the effects of FUV pho- 
tons on the chemistry and structure of the neutral ISM in galaxies, 
as well as diagnosing the conditions within star forming regions. 
PDRs are responsible for most of the infrared radiation from galax- 
ies. The FUV photons usually arise from massive stars creating H II 
regions but sometimes from active galactic nuclei (AGN), which 
produce strong ultraviolet (UV) and X-ray emission. 

Numerical models of PDRs have been around for about 30 
years or so and they have now evolved into complex computer 
codes accounting for a large number o f physical and chemical ef- 
fects (see review by Rollig et al. 2007, hereafter "R07"). In par- 
ticular, the past decade has seen a proliferation of codes capa- 



ble of treating the chemistry and thermal balance within PDRs, 
each developed with distinct interests in mind. Some codes aim 
to encapsulate the detailed microphysics that describe the chem- 
ical an d thermal processes at work in the gas and on the grains 
(e.g. iLe Petit et alj |2006| . 120091 ; iFerland et al.1 [l998l lAbel et al.l 
l2005l 2008), while others focus on treating the gas-grain chem- 
istry or other specific processes in detail whilst approximating 
others in order to explore large regions of p a rameter space (e. 
IWolfire et all 120081 : iHollenbach etai] |2009| ; iRollig et al.l [200 
Some have been developed to model specific source structures 
with either one-plus-one-dimensional or fully two-dimensional 
geometries, including disks and outflow c avities around proto- 
stars (e.g., iKamp & van Zadelhofj 1200 ll : iBruderer et all (20091: 
IWoitke, Kama & Thill2009F >7 further departures from simple one- 
dimensional slab or spherical geometries have been pursued in 
models that treat PDRs as ensembles of discrete clumps, described 
by size and mass dist ribution functions, embedded within an inter- 
clump medium (e.g., ICubick et alj2008l ; lKramer et al.ll2008h . 
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Efforts have also been made in model ing detailed micro- 
physics in dynamically evolved simulations. iGlover & Mac Low] 
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d2007dlbh and lDobbs etal] d2008h have included the formation of 
H2 in three-dimensional simulations of clo ud formation. Improve- 
ments of these methods have been made bv lGlover et al.l ( l2010h to 
additionally model the CO formation thus paving the way towards 
the comparison between simulations and observations. 

Observations of atomic fine structure and molecular lines from 
PDRs have improved with the advent of infrared (e.g. ISO and 
Herschel) and submillimeter (e.g. JCMT and IRAM) telescopes, 
while models have been benefiting from increasingly accurate lab- 
oratory and theoretical data. Most models feature plane-parallel ge- 
ometry, illuminated on one or both sides. This simplifies the ra- 
diative transfer problem because the illumination comes from one 
side only and hence only one line of sight needs t o be taken into 
consideration jFlannerv. Roberge, & Rvbickilll980l R07 and ref- 
erences therein). In addition, some models use a spherical geome- 
try where an isotropic FUV irradiation is taken into consideration. 
R07 provide a detailed account of the differences between plane- 
parallel and spherical models and underline the many assumptions 
and approximations implicit in both geometries, in particular when 
it comes to the treatment of the attenuation due to dust. 

Neither of these two approaches can however deal with more 
complex geometrical issues such as, for example, dumpiness in- 
side the clouds, multiple sources of radiation and non-spherical ge- 
ometries of H II regions or galaxies. For extragalactic sources in 
particular, often unresolved at far-IR wavelengths, accounting for 
the total emission in fine structure and molecular lines, as well as 
from dust, requires the modeling of ensembles of star formation 
complexes where the geometrical issues raised above become im- 
portant. Their modeling is best achieved by simultaneously mod- 
eling the observed spectra and structures of H II and PDR com- 
plexes. Hence three-dimensional integrated photoionization and 
PDR codes are essential for interpreting the wealth of data avail- 
able for star-forming galaxies. However, no code currently offers 
a complete three-dimensional treatment of both ionized and PDR 
regimes (e.g. with realistic radiation fields and geometries and mul- 
tiple exciting sources). Such a self-consistent code is particularly 
important for the modeling of external galaxies where the avail- 
able angular resolution is not high enough to disentangle the dif- 
ferent gas components. While three-dimensio nal gas and dust pho- 
toionization codes exist, such as M OCASSIN JErcolano et al . 2003; 
lErcolano. Barlow, & Storevll2005h . a fully three-dimensional code 
for PDRs that can handle the gas-grain chemistry as well as the 
thermal balance is still lacking. 

In this pape r we present a development of the UCL_PDR code 
( iBell et al]l2006l) to treat three-dimensional structures of arbitrary 
density distribution. The UCLJPDR code is an already benchmarked 
(R07) one-dimensional time- and depth- dependent gas-grain PDR 
code which includes time-varying density and radiation profiles. 
Our new code, 3D-PDF0, adopts the same features in modeling 
the chemistry as UCL.PDR does. It is a starting point towards the 
implementation of an integrated code which will treat dust, pho- 
toionized gas and PDRs together in fully three-dimensional com- 
putational domains. The integrated code aims to couple 3D-PDR 
and MOCAS SIN, with the latter treating the propagation and attenu- 
ation of the UV radiation field as realistically as possible, including 
a detailed spectral energy distribution (SED) profile. 

The paper is organized as follows. In Section[2]we discuss the 
numerical treatment, giving an overview of the code, our ray trac- 
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ing scheme, treatment of the escape probability method, treatment 
of gas cooling and gas heating, the thermal balance and conver- 
gence criteria, as well as the approximations and assumptions we 
made. In Section [3] we present angular and spatial resolution tests 
for the requirements of our ray tracing scheme, and we benchmark 
our code against various one-dimensional codes discussed by R07. 
In Section [4] we show three examples to demonstrate the capabili- 
ties of our code in simulating one- or two- component UV fields in 
uniform or arbitrary density distributions. We discuss our conclu- 
sions in Section[5] 



2 NUMERICAL TREATMENT 
2.1 Overview 

The 3D-PDR code uses the chemical model fe atures of the full y 
benchmarked one-dimensional code UCL_PDR llBelletalJ f2006). 
It solves the chemistry and the thermal balance self-consistently 
within a given three-dimensional cloud of arbitrary density distri- 
bution. We note that 3D-PDR has the ability to solve any one-, 
two- or three- dimensional structure, however in this paper we will 
present calculations from one- and three- dimensional PDRs only. 
The code uses a ray-tracing scheme based on the HEALPix pack- 
age (see Q2.2\ to calculate the total column densities and thus to 
evaluate the attenuation of the far-ultraviolet (FUV) radiation into 
the region (see $2.31 , and the propagation of the FIR/submm line 
emission out of the region. An iterative cycle is used to calculate the 
cooling rates (see $2.5\ using a three-dimensional escape probabil- 
ity method (see $2A\ , and heating rates (see Q2.61 . At each element 
within the cloud, it performs a depth- and time- dependent calcula- 
tion of the abundances for a given chemical network (see $2.71 to 
obtain the column densities associated with each individual species. 
The iteration cycle terminates when the PDR has obtained thermo- 
dynamical equilibrium, in which the thermal balance criterion is 
satisfied (see t]2.8t i.e. the heating and cooling rates are equal to 
within a user-defined tolerance parameter. 

In $2.9\ we present the approximations and assumptions we 
made for the three-dimensional treatment of PDRs and in Appendix 
lAl we present a flowchart of the computational scheme used in 
3D-PDR. 



2.2 Ray tracing 

The three-dimensional ray-tracing scheme we use for calculat- 
ing the column densities and the attenuation of t he FUV ra- 
diatio n field is based on the HEALPix algorithm (Gor ski et al.l 
2005). HEALPix has been u s ed in the past for similar purposes 



(i.e. lAbel & Wandeltj I2OO2I ; lAlvarez. Bromm. & Shapiro! 
Abel, Wise, & Bryanl 120071 iKrumholz. Stone. & Gardined 



2006; 



2007; 



Bisbas et al.ll2009l ; IClark. Glover. & Klessenll2012h . It creates a set 
of Mi = 12 x 4 e pixels uniformly distributed over a unit celes- 
tial sphere, where I is the level of refinement. Each of these pixels 
represents the end of a vector (hereafter 'ray') emanating from the 
centre of a Cartesian co-ordinate system. A pixel defines the centre 
of an approximately square element of solid angle fie — 4-k/Mi. 

Consider a cloud with an arbitrary density distribution consist- 
ing of Aricm elements. Let p(x,y,z) be a random element from 
which the HEALPix rays are emanated all over the computational 
domain. Thus the cloud will be divided into sub-volumes, the ele- 
ments of which belong to different rays of solid angle Qi . For eval- 
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Figure 1. This figure illustrates how evaluation points are created in 
3D-PDR. Gray-filled circles are the elements of the cloud. The white cir- 
cle on the left represents the element p from which a HEALPix ray (solid 
line) emanates. Dashed lines show the extent of solid angle f2£. Black dots 
are the evaluation points. Dot-dashed lines show the extent of the search 
cone which has as vertex the fc-th evaluation point and apex angle 20 CI .; t 
(0 < 9 cr it ^ vr/2 is user-defined). The projection of an element p; on 
the HEALPix ray will be taken if (j>k i _ 1 = Pi^i—iq ^ $crit. where q is 
the HEALPix pixel, creating a new evaluation point. Every new evaluation 
point defines the vertex of the new search cone which however keeps the 
same apex angle in the sense that the cone 'moves' in parallel as we walk 
along the HEALPix ray. 



uating the integrations along each of these rays we create a discrete 
set of points which we dub 'evaluation points'. 

The evaluation points are created by projecting the elements of 
the cloud which are closest to the line of sight of a specific ray (see 
Fig[TJ. Similar schemes for creat ing evaluation points have been 
used also by other workers (i .e. iDale. Ercolano. & Clarke! [2007; 
iKessel-Devnet & Burkert|2000l) . The steps we follow are described 
below. 

(i) We first sort all A/eiem elements with increasing distance 
from the element p using a heapsort algorithm. 

(ii) For a random element p' and using the ang2pix_nest 
HEALPix subroutine, we find the pixel q in whose solid angle the 
element p' belongs to. 

(iii) If the angle <f> p = p'pq is <j> p =SC # cr it, where < # cr it ^ 
7r/2 is a user defined critical search angle, we take the projection of 
p' onto the respective ray. This projection is the evaluation point ki 
and we assign it with identical properties to those of the projected 
element p . The search angle 9 CI n defines a search cone with apex 
angle 2# cr it and solid angle f2 C rit = 27r ( 1 — cos 6> cr it ) whose vertex 
is, in this case, the element p. 

(iv) The evaluation point ki defines the next angle, 0^ , between 
another element, p" , and the pixel, q, i.e. cj)^ = p"kiq. ki also 
defines the new vertex of the search cone, which however keeps the 
same, 2# cr it, apex angle value in the sense that the cone 'moves' in 
parallel as we walk along the ray. 

We repeat the above steps until we reach the end of the ray, for ev- 
ery HEALPix ray; and for all A/'eiem elements considering that each 
one of them is a HEALPix source. We store all of these hierarchies 
in memory. 

For a given spatial function f(r) (e.g. number density of 
some chemical species), we perform integrations along each ray 



by adopting the trapezoidal rule: 

f L a u "v^ /( r fe-i) + f( r *) i i m\ 

/ /(r)<2r~ 2^ k \ r *>~ rk - 1 \- (1) 

Jo fc=i 

Here, L is the length of the ray, A4val is the total number of evalu- 
ation points along this ray, and is the distance of the evaluation 
point k from the element p which defines the HEALPix source. 
The distance |rjt — rk-\\ is equivalent to the spatial distance be- 
tween two evaluation points of the same ray which in turn defines 
the integration step dubbed as an 'adaptive step'. 

The advantage of this ray-tracing scheme is that it can be ap- 
plied directly to both grid-based and Smoothed Particle Hydrody- 
namics data without the necessity of implementing further modifi- 
cations. 



2.3 Treatment of the ultraviolet radiation field 

For the purposes of this paper and for the current verion of 
3D-PDR, we simplify the treatment of the UV field; we neglect the 
contribution of t he diffusive radia tion by invoking the on-the-spot 
approximation dOsterbrocklll974l) . We do this in order to explore 
the effects introduced when one moves to a three-dimensional treat- 
ment of PDRs. As mentioned in the Introduction, o ur future plans 
include the coupling of 3D-PDR with M OCASSIN dErcolano et al.l 
l2003l;lErcolano, Barlow. & StorevteOoH) in a single integrated code 
in order to include a realistic treatment of a three-dimensional UV 
radiation field, and that is including the diffusive radiation compo- 
nent. In the current version of 3D-PDR the user is able to choose 
between three types of UV radiation field: plane-parallel (UNI), 
radial-sampling (ISO), or a field emitted spherically symmetrically 
by a point source (PNT); or any combination between these. 

For a random element p(x, y, z), t he strength of the incident 
radiation, x(z>). measured in units of the lDraineR l978) interstellar 
radiation field, is calculated using the equation 

Xip) = f o 4 \o(9,<t>)e-™ Av{M) § 

~ ^Exo(q) C - Tuv ^ (q \ ( 2 ) 

q=0 

where xo(#, 4>) is the magnitude of the unattenuated field strength 
(Draines) at the surface of the cloud in direction 6 and <j>, and 
i"uv = 3.02 is a dimensionless factor converting the visual extinc- 
tion to UV attenuation. The integration over solid angle ( J dQ.) is 
approximated using a summation over all Ne HEALPix rays. The 
term Av is the visual extinction defined as 

Av(<l) = A Vo I n H dr 
Jo 

* A Vo J2 n ^ k -V + nH{k) Ar, (3) 

k=l 

where Av a ~ 6.29 x 10 - 22 mags cm 2 . The integration (f n H dr) 
corresponds to the column of the H-nucleus number density n H 
along the ray of length L. This integration is approximated using a 
summation over the M CV ai evaluation points of the q HEALPix ray 
and where Ar = |r* — rk-i\ is the adaptive step as discussed in 
Eqn.lfTJl. 



4 T. G. Bisbas, T. A. Bell, S. Viti, J. Yates and M. J. Barlow 



2.4 Escape probability 

Suppose that the random element, p(x, y, z), of the cloud is con- 
sidered as a source of radiation. Assuming statistical equilibrium, 
the level populations and radiation field are related by 



* (P) E Ri J ( p ) = E n i w ' 



(4) 



where the summation is over the total number of levels included 
and m (p) , rij (p) are the populations of levels i, j respectively. The 
left hand side describes emission and the right hand side absorption. 



The term Rij (p) is 



Aij + Bij(Jij(p)) + Cij(p), i > j 
Bij(Jij(p)) +Cij(p), i < j 



(5) 



where Aij and Bij are the Einstein coefficients, Cy(p) is the col- 
lisional rate for excitation (i < j) and de-excitation [i > j), 
and {Jij(p)) is the mean integrated intensity received at p from 
all solid angles dfi. In our models, we adopt the Large Velocity 
Gradient (LVG) or escape probability formal ism ( Sobolev 1960; 
Castodll970k lde long. Dalgarno. & Chulll975l : iPoelman & Spaans 
2005) to describe the mean radiation field { Jij (p)) as: 



(Jy (p)> = [1 - Pij(p)]Sij(p) + Pij(p)B(v k 



(6) 



The term Sij (p) is the source function due to transitions between 
levels i, j and is given by 



Sij(p) 



2hvf 



n t (p)gj 



(7) 

gj are the statistical 



where u%j is the photon frequency, and 
weights of rij , rii respectively. 

The term B(yij) is the total background radiation valid at 
FIR and submm wavelengths, including Cosmic Microwave Back- 
ground blackbody emission at T C mbr = 2.7 K and dust emission 
approximated as a modified blackbody at Td ust . 

The term /9y (p) describes the probability that a photon of fre- 
quency Vij escapes from the element p without interacting with the 
rest of the cloud. In the present work we consider wavelengths in 
the FIR and submm range and we therefore neglect the absorption 
due to dust, considering only the line absorption component. We 
adopt the analytical expression for the esca pe probability fiij de- 
veloped bv lde Jong. Dalgarno. & Chul dl975t) : 

Pa 



r 471 dn 


"1 - e" TL " 


L 4^ 





(8) 



where r L = T%j (p, q) is the line optical depth at the element p 
along the direction q, given by the expression 



Tij(p,q) 



8^ 3 



Au(p) 



nj(p)gi 



dr , 



(9) 



where the integration is performed between the positions fj and r 2 
which define the direction q, and Au(p) is the root mean square of 
the thermal and turbulent velocities (see below). 

We approximate the integral over all space of Eqn.([8j using 
the HEALPix ray scheme described in ij2.2l as 



/ ~~a — = T7~ / HEALPix 



(10) 



Thus for the element p, the numerical expression for the escape 



probability is 



1 

M 



Mi 

E 

q=0 



j (p.q) 



T«(p,q) 



(11) 



where q represents each individual HEALPix ray. This equation 
provides the total escape probability corresponding to the summa- 
tion of all individual escape probabilities per HEALPix direction 
averaged over the total number of these directions. The numerical 
expression of Eqn.l[9]l for the line optical depth, Tij (p, q), is 



Tij(p>q) 



8nufjAu(p) 



A/eval 

E 



[rij(k-i)+nj(k)}gi 



2g 3 

m{k — 1) + nj(fe) 



Ar, (12) 



where the summation is over all A/" cva i evaluation points along the 
q ray, Ar = \ru — rjt_i| is the adaptive step, and Aii(p) is given 
by 



Au(p) = ^ 



8k B T{p) 
■KTUn 



(13) 



where fc B is the Boltzmann constant, m H is the proton mass, T(p) 
is the gas temperature of element p and ii tbrb is the user-defined 
turbulent velocity. 

The LVG or escape probability method described here is used 
as an approximation to describe the mean radiation field intensity. 
Our future plans include the treatment of an exact line transfer from 
one region to another inside the PDR. 



2.5 Treatment of gas cooling 

The gas in molecular clouds is cooled primarily by the collisional 
excitation and subsequent emission of a number of key atomic and 
molecular species. Within PDRs this is usually dominated by emis- 
sion in the [C II], [O I] and [C I] fine-structure lines and in the ro- 
tational transitions of CO. The cooling rates due to emission from 
these species are determined by the 3D-PDR code at every element 
within the cloud by calculating the emissivity (in erg s" 1 cm ) 
of each transition, having solved for the non-LTE excitation and 
radiative transfer under the LVG assumption, as described in i]2,4l 
The total radiative cooling rate at each position within the cloud is 
then the sum of the emissivities of all radiative transitions consid- 
ered. We include transitions between the ground-state fine structure 
levels of O ( 3 P 2 , 3 A, 3 P ), C ( 3 P , 3 A, 3 P 2 ) and C+ ( 2 P 1/2 , 
2 an d 1 1 rotational levels of CO (note that up to 40 rotational 
levels can be treated in the code, but we limit the number here to 
increase computational speed for these initial models). Collisional 
excitation rates are t aken from the Leid en Atomic and Molecular 
Database (LAMDA; Ischoier et alj|20051) for all available collision 
partners, namely H, He, H 2 , H + and e~ for O and C; H, H 2 and 
e~ for C + ; and H, He and H 2 for CO. The collisional rates at the 
required gas temperature are determined by linear interpolation be- 
tween the fixed temperature values specified within these data files. 

Deeper within molecular clouds, other molecular species, in- 
cluding the isotopologues of CO, OH and H 2 0, become important 
coolants, but despite being included in the UCL.PDR code we ne- 
glect their contribution within the 3D-PDR code since we are con- 
cerned with the thermal balance near the surfaces of these clouds, 
where gas temperatures show the greatest variations. At high vi- 
sual extinctions, where these other molecular coolants become im- 
portant, the gas and dust temperatures generally tend to a rather 
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constant 8-15 K (e.g. Fig. 12 in R07) in the absence of embedded 
heating sources and the appropriate dark cloud chemistry asserts 
itself, with little sensitivity to these small temperature variations. 

At high densities (> 10 6 cm -3 ), collisions with dust grains 
can also efficiently cool or heat the gas, depending on the tem- 
perature difference between the gas and the grains. This mech- 
anism is also accounted for in the 3D-P DR code, following the 
treatment of Burke & Hollenbach dl983h and the accommoda- 
tion fitting formula of iGroenewe gen ( 1994), assuming a standard 
iMathis. Rumpl & Nordsieckl (MRN: [l977l) grain size distribution. 

Cooling due to ro-vibrational emission of H2 can also play 
a minor role in regulating the gas temperature close to the cloud 
surface, but we do not include this in the current version of the 
3D-PDR code. At relatively high temperatures (> 5000 K) neutral 
atomic gas can be efficiently cooled by excitation of the metastable 
1 D levels of atomic carbon and oxygen. These processes are im- 
plemented in the code, but are disabled for the present study. 



2.6 Treatment of gas heating 

The mechanisms that contribute to the total heating of the gas and 
their treatment in the 3D-PDR code are identical to those in the 
UCL_PDR code, described in detail bv lBelletal.ld2006l) . Near the 
PDR surface, the dominant gas heating mechanism is the pho- 
toelectric ejection of electrons from small dust grains and poly- 
cyclic aromatic hydroc arbons (PAHs). We adopt the treatment of 
iBakes & Tielensl dl994l) . where the total heating rate (in units of 
erg cm -3 s -1 ) is given by 



10 eG n H , 



(14) 



where 



4.87 x 10" 



[1 + 4 x lO-^GoTVz/ne) - 73 ] 
3.65 x 10 -2 (T/10 4 )"- 7 



+ 



(15) 



[1 + 2 x 10- 4 (G TV2/n e )] 

and G p is the local FUV flux expressed in units of the |Habing 
dl968l) field (which is related to the scaling factor for the Draine 
field by Go = l-7x)> is the total H-nucleus number density 
(cm -3 ), T is the gas temperature (K) and n e is the electron number 
density (cm -3 ). This heating rate is countered at high gas tempera- 
tures by cooling due to recombination of electrons with grains and 
PAHs, and we sim ilarly adopt the analytical expression derived by 
IBakes & Tielensl J 19941) for this cooling rate: 

A REC = 3.49 x 10 -30 T ' 944 (GoT 1/2 K)Vn H (16) 
P = 0.735/T ' 068 . (17) 

We assume the same standard MRN grain size distribution adopted 
by those authors in deriving these rates. We note that Eqn.dl4t de- 
scribing the universal grain heating rate should be used for low in- 
tensities of radiation field such as those in th e simulations presented 
here. In future updates of 3D -PDR the IWeingartner & Draind 
heating rate will be used which gives better approximation 
at higher intensities of radiation fields. 

The collisional de-excitation of vibrationally excited H2 fol- 
lowing FUV pumping can also be an important heating mechanism 
in dense gas near the cloud surface. We assume a single excited 
pseudovibrational level of H2, denoted H2, to effectively account 
for the full distribution of H2 molecules in vibrationally excited 



levels and we describe our treatment of its formation and destruc- 
tion in the next section. We adopt the associated heating rate fo r 
collisional de-excitation of H2 from lTielens & Hollenbach! dl985h : 

Th* = [n(H) 7 » H + n(H2)7* H oXH2),E*, (18) 

where n(H), n(H2) and n(H|) are the number densities of H, H2 
and H2, respectively (cm -3 ), E* is the energy of the single excited 
pseudovibrational level (2.6 eV), and 7^ and 7„ 2 are the colli- 
sional de-excitation rate coefficients (in units of cm s - ) from the 
excited to the ground vibrational level for H and H2, given by 



H 

7»o 

Y*o 



10 -12 T l/2 e -1000/T 

1.4 x io -12 r 1/2 e -18100/(T+1200) . 



(19) 
(20) 



We note that more complicated techniques for the H2 treatment 
have been implemented and modelled showing differences in the 
heating rate of up to ~ 1 dex (see R07). 

In addition, photoi onization of n eutral carbon liberates about 
1 eV per photoelectron dBlack|[l987l . the rate of carbon photoion- 
ization is described in the next section). The internal energy of 
newly-formed H2 as it leaves the grain surface can also make a 
non-negligible contribution to the gas heating near the PDR sur- 
face; following Black & Dalgarno (1976), we assume that the 4.48 
eV of internal energy is distributed roughly equally between trans- 
lation, vibration and rotation, so that 1.5 eV will go into kinetic 
energy per H2 molec ule formed. D eeper within the cloud, cosmic 
rays and turbulence telackj |T987h do most of the heating of the 
gas , with exothermic reac t ions a lso playing a minor role. Follow- 
ing lTielens"& Hollenbach ( 1985), we assume a cosmic ray heating 
rate of 



r CR = 1.5 x 10 -11 Cn(H 2 ), 



(21) 



where ( is the cosmic ray ionization rate per H2 molecule. For the 
rate of g as hea ting due to dissipation of supersonic turbulence, we 
assume dBlackll 19871) 

Tturb = 3.5 X 10 ^turb/^turb^h, (22) 

where d T urd is the turbulent velocity (km s -1 ) and £tur.b is the 
turbulent scale length (assumed here to be 5 pc). 



2.7 Model chemistry 

The code determines the relative abundances of a limited number of 
atomic and molecular species at each cloud element in the model 
by solving the time-dependent chemistry of a self-contained net- 
work of formation and destruction reactions. The chemical net- 
work i s a subset of the m ost recent UMIST database of reaction 
rates dWoodalletaT]|2007l) . consisting of 320 reactions between 
33 species (including electrons), and includes photoionization and 
photodissociation reactions in addition to the standard gas-phase 
chemistry. We make use of the XDELOAD tool (kindly provided by 
L. Nejad) to construct the set of ODEs describing the formation 
and destruction of each species and the associated Jacobian matrix 
that together are used to compute the chemical abundances in the 
3D -PDR code. We order the species in the ODE system according 
to their number of formation reactions, which has been sh own to 
significantly speed up the computation of the abundances dNeiaj 
2005). In this paper we restrict our models to produce steady-state 
abundances by assuming a chemical evolution time of 100 Myr, 
sufficient for all reactions to reach equilibrium. However, the code 
is capable of following the full time-dependent evolution of the 
chemistry within a cloud, making it a powerful tool for modeling 
dynamically evolving structures. 
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In our full chemical network, we include reactions involv- 
ing vibrationally excited molecular hydrogen, whose internal en- 
ergy can provide a means to overcome the activation barrier of 
certain neutral-neutral and ion-molecule reactions, thus consider- 
ably enhancing the abundances of some species, such as CH + . 
Significant abundances of molecular hydrogen in vibrationally ex- 
cited levels can be maintained in PDRs due to the FUV pump- 
ing of H2 to electronically excited states, followed by its decay 
to populate the vibrational levels of the ground electronic state. 
In our models, we adopt a simplified trea tment of this mecha- 
nism following riTelens & Hollenbachl dl985h and others in assum- 
ing that all vibrationally excited H2 is in a single characteristic vi- 
brational level, with effective spontaneous emission and collisional 
de-excitation rates, that approximates the full distribution amongst 
all vibrational levels. Molecular hydrogen in this pseudo-level is 
labell ed H2 and the reactions tha t form and destroy it are taken 
from lTielens & H ollenbach ( 1985), with updated rates for some re- 
actions Ttakenfix)rri|Ag^nde^ray d2010h . 

We also include reactions involving neutral and singly ionized 
polycyclic aromatic hydroc arbons ( PAHs ) , adop ting the rate coef- 
ficients proposed by Wol fire et al.l J2003l 120081) . though we omit 
these reactions for the reduced chemistry of the benchmark com- 
parisons and test applications described in this paper. More detailed 
models to be presented in forthcoming papers will include the full 
reaction network, including the role of PAHs and vibrationally ex- 
cited H2 . 

The formation of molecular hydrogen on grain surfaces is a 
key reaction in determining the properties of PDRs since, together 
with destruction through photodissociation (described below), it 
governs the transition from atomic H to H2 within the PDR and 
therefore plays a critical role in the chemistry (many species form 
through reactions with H2) and thermal balance (through heating 
processes such as those involving vibrationally excited H2, and as 
a collision partner for coolant species and a coolant in its own 
right). We h ave implemented the treatment of ICazaux & Tielensl 
(2002, 120041) who have shown that the process can be decribed us- 
ing a rate equation formalism that adequately accounts for both 
chemisorbed and physisorbed hydrogen atoms reacting on grain 
surfaces. We use their standard values for the properties of both 
graphitic and silicate grains to determine the formation efficiency, 
and t he expression for the sticking coefficient of H atoms on grains 
from|Hollenbach & M cKeel ( fl979l) . We describe the grains by their 
global population properties, assuming a standard MRN size dis- 
tribution and graphite/silicate composition. We do not consider the 
level-specific distribution of newly-formed H2 leaving the grains, 
instead treating it as one species. For the benchmark tests described 
in E13.3I we have use a simplified formation rate (as described in 
that section) in order to better match the benchmark model specifi- 
cations. 

With the exception of the reaction rates for the photodisso- 
ciation of H2 and CO and the photoionization of carbon, all pho- 
toreaction rates are calculated using the standard UMIST treatment 
adapted to our HEALPix based ray-tracing scheme, where the total 
rate at a given position in the cloud is the sum of the rates deter- 
mined along each HEALPix ray to the PDR surface. The rate (s _1 ) 
for a photoreaction i along a particular ray q is then given by 

ft(q) = a ao (q)e- 7 ^ (q) , (23) 

where a; is the unattenuated photoreaction rate (s _1 ), evaluated for 
isotro pic illuminati on by the standard Draine interstellar radiation 
field Jpraind[l978h . Xo(q) is the FUV flux incident on the PDR 
surface at the element intersected by ray q and specified by a scaled 



equivalent of the Draine field, Av(<l) is the total visual extinction 
along that ray to the surface, and 7, is a scaling factor that relates 
the attenuation in the visible to that in the FUV. 

The photodissociation rates of H2 and CO depend sensitively 
on the column densities of these species along the direction of 
the incident FUV radiation, since their absorption of FUV pho- 
tons leads to significant shielding against photodissociation. This 
(self-)shielding is therefore treated explici tly in the code, using 
the results of the detailed c alculations of iLee et ID lfl996h and 
Ivan Dishoeck & Black! J 19881) for H2 and CO, respectively. Those 
authors ran full radiative transfer models accounting for both self- 
shielding and line overlap of H, H2 and CO lines in order to de- 
termine the degree of shielding produced under a range of cloud 
conditions. They found that these detailed processes could be 
well described by shielding functions that depend only on the to- 
tal H2 and CO column densities to the PDR surface. We have 
therefore adopted these treatments and use the tabulated shield- 
ing functions provided in their papers. We do not explicitly ac- 
count for state-specific photodissociation of H2 and CO, since 
the rates and self-shielding fa ctors listed by ILee et al] i ll 9961) and 
Ivan Dishoeck & Black! J 1988t) are strictly valid for the global pop- 
ulations of these molecules, with the exception of our inclusion of 
vibrationally excited H2 in a single excited pseudo-lev el, for which 
we ad opt the effective photodissociation rate given by Rol lig et"al] 
(2006). We also neglect detailed treatment of ro-vibrational cas- 
cades following electronic excitation by the UV photons, since such 
a treatment would dramatically increase the computational time 
without significantly altering the resulting abundances. In addition, 
we account for the shielding of neutral carbon again st photoion- 
ization using the treatment o fKamp&Bertoldi(2000). In all three 
cases, the column densities of H2, CO and C needed to calculate 
the shielding factors at a given point in the cloud are determined 
along each HEALPix ray to the PDR surface. 

The calculation of the H and H2 abundances at the elements 
near the cloud surface depends critically on the shielding provided 
by H2 against its own photodissociation, which, taken together with 
H2 formation on grains, represents the main formation/destruction 
cycle for molecular hydrogen in the UV-illuminated gas. The 
amount of self-shielding along a given ray to the cloud surface is 
itself sensitive to the total column density of H2 and therefore re- 
quires that the H2 abundance be known at each evaluation point 
along the ray. A similar relation links the photodissociation of CO 
to its column density along each ray. 

There exists, then, an interdependence between the photodis- 
sociation rates and the abundances of all elements near the cloud 
surface, requiring that the abundances be calculated and the result- 
ing column densities updated a number of times before correct val- 
ues can be obtained for both. We therefore perform a chemistry 
iteration each time that the gas temperatures are changed, in which 
the reaction rates and shielding factors are determined at the new 
temperature and for the current column densities, new abundances 
are calculated, the column densities are updated based on the new 
abundances, and the process is repeated for Ichem iterations. After 
numerous tests for convergence, we find that between 5-10 chemi- 
cal iterations are needed at the start of the code in order to correctly 
determine the abundances of H and H2, and C + , C and CO, near 
the surface. Following this first determination of the chemistry at 
the initial "guess" temperature (see i j2.8l >. we find that subsequent 
changes to the gas temperature require only 1 or 2 iterations of the 
chemistry to reach convergence. For the models presented in this 
paper, we have performed Ichem = 8 chemical iterations at the 
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start of the code and then 3 chemical iterations after each change to 
the gas temperature. 

The capability exists within the 3D-PDR code to include the 
role of grains in the chemistry, including freeze-out of gaseous 
species onto grains, surface reactions and the release of grain man- 
tle species back into the gas phase by means of evaporation, pho- 
todesorption or desorption due to cosmic ray heating. However, in- 
clusion of the full gas-grain chemistry can dramatically increase 
the computational time needed to follow the evolving abundances, 
so reactions involving dust grains are currently omitted from the 
models in order to increase the speed of the code and to be able to 
compare with other PDR models that excluded grain chemistry for 
the benchmarking tests of R07. 



2.9 Approximations and assumptions 

In order to make tractable the problem of simultaneously calculat- 
ing the chemistry, thermal balance and radiative transfer within a 
three-dimensional cloud, we have necessarily made a number of 
simplifying assumptions and approximations in the 3D-PDR code. 
We list here the most important of them: 

(i) Assumption: the level populations and resulting emissivities 
change rapidly with respect to the changes in abundance due to the 
chemistry. 

(ii) Approximation: the radiative transfer of the FUV radiation 
into the model cloud is treated by considering only the attenuation 
by dust and neglecting the line scattering/di ffusive term s by invok- 
ing the on-the-spot approximation ( Oster brockl 1 19741) . The grain 
temperature is not affected by the UV radiation field. 



2.8 Thermal balance and convergence criteria 

3D-PDR starts the calculations by setting up an initial "guess" tem- 
perature profile which is used in order to begin the iterative process 
for level populations and thermal balance convergence. The im- 
plicit assumption of an initially uniform temperature profile over 
the entire PDR of arbitrary density distribution may lead to un- 
physical values of level populations in certain parts of the PDR 
causing the iterative process to fail (i.e. level populations cannot 
converge). To avoid this we ran several one-dimensional calcula- 
tions of uniform density PDRs interacting with several FUV field 
strengths. The density of PDRs in these calculations spans 10 2 ^ 
n (cm -3 ) ^ 10 6 and the field strength spans 1 ^ \ (Draines) ^ 
10 6 . We find that the equation 

T gucss = 10[1 + (100 X ) 1/3 ] (24) 

provides an acceptable estimate for an initial "guess" temperature 
profile (first iteration over thermal balance) in order to begin the 
overall iterative process. In this equation, T gucss (K) is the temper- 
ature and x (Draines) is the attenuated FUV field strength. 

Using the temperature profile of Eqn.J24b the code begins the 
iteration process in order to obtain level population convergence. 
We assume that the level populations have converged when the 
change in any given population between two consecutive itera- 
tions is less than a user-defined tolerance parameter (in this paper 
iT crr < 1%). We then calculate the total cooling (A) and heating (T) 
rates. By comparing these rates we assign new gas temperatures i.e. 
if A > r a lower temperature than T guoss is required and if A < T 
a higher temperature than T guoss is required. The technique we use 
is the following. For a cloud element p, if its temperature changes 
monotonically from that given by EqnJ24t. the new temperature 
used in the next iteration over thermal balance will differ by 30% 
from its current value. If the change is not monotonic, then a binary 
chop routine is performed between the current temperature and the 
one obtained in the previous iteration on thermal balance. 

Once the new kinetic temperatures have been calculated and 
updated, we start iterating again over the level populations. We re- 
peat this process (i.e. assigning new temperatures, iterating over 
level populations, etc.) until we reach convergence over thermal 
balance. We assume we have obtained thermal balance conver- 
gence when the heating and cooling rates are equal to within some 
user-defined error tolerance (in this paper a CIItT sC 0.5%) or 
when the change in gas temperature between iterations is negli- 
gible, i.e. smaller than a user-defined, Tdjg, value (in this paper 
Idiff sj 0.01 K). 



3 BENCHMARKING 

3.1 Angular resolution and search angle 

In this section we explore the values of the level, £, of angular res- 
olution and the search angle, 9 C nt, for which 3D -PDR evaluates 
integrations at reasonable accuracy along the HEALPix rays. To do 
this we perform a test to measure the integration accuracy of the 
code in the evaluation of column density. We consider a spherically 
symmetric cloud, the centre of which defines the centre of the co- 
ordinate system. The radius of the cloud is R — 1 pc and its spatial 
H-nucleus number density profile is nn(r) = noe~ r ^ R , where r is 
the radial distance from the centre and no = 100 cm -3 . Therefore, 
the column density, N, in any direction as seen from the centre is 

f R 

N = / n H (r)dr ~ 1.95 x 10 20 cm~ 2 . (25) 
Jo 

For the construction of the cloud we use A4lem = 10 5 ele- 
ments uniformly distributed. We also assume that the hierarchy of 
rays is emanated from the element po which is positioned at the 
centre of the sphere. 

Since the escape probability function /2y (p) of the element p 
(Eqn. lilt is the average value of the escape probabilities over all 
HEALPix rays, it is important to know the integration accuracy of 
the present method in calculating a function averaged over all rays. 
For this reason, we define as 




the average value of column densities over all HEALPix rays, 
where N q is the column density of the q ray calculated as discussed 
in Eqn.([3j. 

Figure [2] shows the error 

a crr = 100^^1% (27) 

between the analytical value, N, and the calculated (Njv^) of the 
column density versus # cr it and for different levels of refinement, 
I. The solid line represents the errors at I = (Ao = 12 rays), the 
dotted line at I = 1 (N\ = 48 rays), the short-dashed at I — 2 
CA/a = 192 rays), and the long-dashed at I = 3 (A/3 = 768 rays). 

Overall we see that a crr < 3% even for # C rit as low as 7r/12, 
so the accuracy is acceptable. We note that by increasing Af e i em , 
the error o CII is decreased. We also observe that by increasing the 




Figure 2. This figure plots the error cr err (Eqn l271 between the analytical, 
N, and the calculated, (Njvv), values of column density versus the critical 
search angle, cr it- The solid line represents the errors at i = (A/o = 12 
rays), the dotted line at £ = 1 (Afi = 48 rays), the short-dashed at i = 2 
CA/jj = 192 rays), and the long-dashed at I = 3 (A/3 = 768 rays). For 
£ = and 

$crit 5 s 7r/6 we achieve reasonable integration accuracy and at 
low computational cost. 



level, £, of angular resolution refinement the error a CII is decreased, 
however at the expense of computational time since more evalua- 
tion points are created. For the tests and applications described here 
we will use £ — and 9 CI n — 0.8 ~ 7r/4. 



3.2 Resolution along a ray 

Here we examine the resolution requirements needed along a 
HEALPix ray in order to establish when our calculations are con- 
verged. We use a one-dimensional cloud of uniform H-nucleus 
number density n H consisting of A4lem elements. Since the density 
is constant, from Eqn.([3]l the length L of the cloud will be given by 
the equation 



L(cm) = A\ 



1.59 x 10 21 (cm" 2 ) 
n H (cm~ 3 ) 



(28) 



where we set Av,m ax = 10. The elements are aligned with two 
opposite HEALPix rays, namel)[f] the rays with ID = 4 and 6 
of I = 0. Consequently, the search angle criterion is eliminated 
and the elements pre-define the evaluation points. Considering that 
these two rays define the x — axis of a Cartesian co-ordinate system, 
we apply a UV field of strength \ from the — x side; that is from 
ray ID = 6. For the rest of HEALPix rays we assign very high 
optical depths, implying that the one-dimensional line represents a 
three-dimensional semi-infinite slab. 

We construct the cloud by creating elements logarithmically 
distributed along the x— axis. We use Ma v elements per Av dex 
and with —5 ^ \og(Av) ^ 1. We run three different tests with 
n H = 10 2 cm" 3 and x = 1 Draine (Tl); n H = 10 3 cm" 3 and 



In 3D-PDR we use the NESTED numbering scheme of rays. 



Figure 3. This figure plots the number Ma v of elements per Ay dex 
versus the relative error <r crr between the A^ tranB (where the H/H2 tran- 
sition occurs) and the respective value for Ma v = 200. We find that for 
Ma v 5* 20 we obtain an error of <r crr < 0.5%. 



X = 10 Draines (T2); n H = 10 cm and x = 10 Draines 
(T3). We vary the number Ma v and we measure the value of visual 
extinction A y, trans at which the H/H2 transition occurs. 

Figure [3] plots Na v versus the relative error a alI between the 
j4v,trans and the respective value for Ma v = 200. In all tests we 
find that we obtain convergence for Ma v as low as 20. In addition, 
we find that the minimum value of visual extinction needed in order 
to resolve the H/H2 transition at this error must satisfy the relation 



3.3 Comparison with the other PDR codes 

In order to assess the reliability of the new code, we have run a se- 
ries of models designed to reproduce the benchmark tests described 
in the R07 comparison study. In that paper, results from a number of 
the most widely used PDR codes were compared for a set of models 
in which the capabilities of all codes were restricted to an agreed 
upon (and much simplified) set of parameters and treatments for 
processes such as the chemistry, UV attenuation and heating rates. 
Table [TJ lists the main physical parameters that varied between the 
four models considered; the reader is referred to the R07 paper for 
full details of the other parameters used in the models. We note that, 
whilst we have adopted the same elemental abundances and treat- 
ment for the attenuation of the UV radiation in our tests with the 
3D-PDR code, a number of differences remain between our code 
and those included in the R07 comparison. In particular, we are us- 
ing a slightly larger and updated chemical network with rates taken 
from the most recent release of the UMIST database, and updated 
collisional excitation rates taken from the LAMDA database. Most 
importantly, we are using a modified form of the standard rate of 
H2 formation on grain surfaces (in units of cm 3 s" 1 ), taken from 
Ide Jond ( fl977h . that includes an additional exponential term to re- 
duce the formation efficiency at high temperatures: 
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3 x iq-^t 1 /^-^ 1000 . 



(29) 



As discussed below, this additional term leads to some minor 
differences between our results and those of the R07 benchmark 
models. Whilst this treatment has been superseded in recent years 
by more advanced formalisms that account for both chemisorp- 
tion and physisorption a nd for mixed grain compositions (see, e.g., 
ICazaux & T ielens 20(3), we have chosen to neglect such treat- 
ments for the simplified models we present here. 

To allow comparison with the R07 results, we have restricted 
the cloud geometry used in the models to that of a one-dimensional 
slab of constant density illuminated by a UV field from only one 
side, as described in §3.2. In addition, we have restricted the es- 
cape of photons from the cloud to a single ray directed towards 
the illuminated surface, thereby emulating a semi-infinite slab. The 
gas temperature is allowed to vary and is determined by solving the 
thermal balance explicitly (see §2), but we assume a fixed dust tem- 
perature of Tdust = 20 K throughout the cloud. In all models, the 
cloud consists of Ma v = 100 elements logarithmically distributed 
per A v dex with -4 ^ \og{A v ) < 1.3 for the n H = 10 3 cm -3 
density case (so jVeicm = 530) and with —6 ^ log(-Av) ^ 1.3 for 
the n H = 10 5 ' 5 cm -3 density case (so7V e i cm = 730). 

The results of the 3D-PDR test models are compared 
to those of the R07 paper in Figs|4] [5] and [6] (for discussion 
on Fig[6] see ^4. It . We use t he workshop res ultfl of t he fol- 
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2005; Shawetal 



J l2005h 
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1985; Kaufman ct ; 
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HTBKW jTielens & HollenbacrJ 
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LEIDEN 

r 



19961: iBensch et alJ 
(Black & van Dishoeckl 



van Dishoeck & Black! 119881 : Ijansen et al.1 Il995h, MEIJERINK 
(Meijerink & Spaans 2005) , MEUDON jLe Petit. Roueff. & H erbst 
20041 ; iLe Petit, Roueff. & Le Bourlotl |2002|; iLe Bourlot et all 
19931) . and STERNBERG dSternberg & Dalgarnol Il989l Il995l; 



2003 



1987 



Boger & Sternberg 2006). 



Figure [4] shows the results for benchmark model V2 (n H = 
10 3 cm -3 ; x = 10 5 Draines), including the gas temperature pro- 
file, number densities of H, H2, C + , C and CO, and emergent 
intensities (surface brightnesses) of the dominant cooling lines. 
As can be seen, the overall agreement is very good, with the 
results from the 3D-PDR code typically falling within the scat- 
ter of results produced by the other codes. In addition to the re- 
sults from the R07 comparison study, we have also obtained re- 
sults for the four benchmark models using the latest version of the 
UCL_PDR code by adopting identical physical parameters and the 
same chemical network as used in the 3D-PDR code. The results 
for model V2 are included in Fig|4| (labelled as UCL_PDR11) and 
show excellent agreement between the one-dimensional and three- 
dimensional versions of the UCL_PDR code. The results are simi- 
larly identical for the other three benchmark models and are there- 
fore not shown in the remaining figures. 

The only notable differences between the 3D-PDR and the 

3 NOTE: In order to perform a comparison of 3D-PDR 
with other available PDR codes we use some of the data of 
models VI, V2, V3, and V4 taken from the workshop site 
I http://www.astro.uni-koeln.de/site/pdr-comparison 1. Due to lack of 
convergence or incompleteness we do not always include every code in our 
plots. 



Table 1. The one-dimensional benchmark models performed for compari- 
son with the results of R07. We refer the reader to that paper for full details 
of the model parameters used. 



Model ID 


ra H (cm 3 ) 


X (Draines) 


VI 


10 3 


10 1 


V2 


10 3 


10 5 


V3 


10 5.5 


10 1 


V4 


10 5.5 


10 s 



R07 results visible in Fig|4]are the much lower H2 abundance at 
the outer cloud edge, which is due to the reduced H2 formation 
efficiency at high temperatures, as discussed above, and the rise 
in neutral carbon abundance (and corresponding drop in C + abun- 
dance) at lower Av in the 3D-PDR model, which is due to the more 
advanced treatment of the carbon photoionization rate that we have 
adopted, including shielding by lines of H2 and C (see i]2.8l >. This 
difference is also reflected in the [C I] local emissivity profile. 

Figure [5] shows selected comparisons of the results for bench- 
mark models V3 and V4. Although we state that the benchmark- 
ing model V3 of R07 should not be considered as a potential PDR 
(significant contrast between the high density and the low radiation 
field which leads to a temperature of ~ 20 K at Av < 0.1 mag), 
we include it in the present work as we are able to compare 
3D-PDR with the other codes even under such extreme conditions. 
While the scatter in the range of results from all the codes is larger, 
the results obtained by our code nevertheless continue to show very 
good agreement with the rest of the codes, the only differences of 
note coming from the different treatments for the H2 formation and 
carbon photoionization rates already discussed. 

Overall, the results of these four benchmark tests demonstrate 
that the 3D-PDR code compares very favourably with the well- 
established PDR codes included in the R07 study when using sim- 
ilarly limited chemical networks and treatments of the various mi- 
crophysical processes. We therefore consider the new code to be 
reliable and in the next section we go on to demonstrate some of 
the more advanced model applications that become possible with 
the fully three-dimensional geometry offered by the 3D-PDR code. 

In addition, we note that 3D— PDR is approximately two times 
faster than UCL.PDR for running one-dimensional models. The pri- 
mary source of this significant speed-up is the usage of pre-defined 
evaluation points which are either user-specified (i.e. in case of one- 
dimensional runs of 3D -PDR) or created automatically due to the 
projection of the elements that make up the cloud (as described 
in i]2.2l i.e. in case of two- or three-dimensional runs). Thus the 
evaluation points act as a fixed grid, in contrast with the adaptive 
grid used in UCL_PDR. However, techniques controlling an adap- 
tive increase of resolution inside PDRs are inevitably necessary for 
3D-PDR in treating complex three-dimensional structures. We plan 
to implement and examine these techniques in a forthcoming paper. 



4 APPLICATIONS 

In this section we present three different applications to demon- 
strate the capabilities of our code in simulating three-dimensional 
cloud structures. We explore i) a uniform-density spherical cloud 
interacting with a plane-parallel external radiation field ( i]4.1b , ii) a 
uniform-density spherical cloud interacting with a two-component 
external radiation field ( £]4.2t . and iii) a cometary globule inter- 
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Figure 4. Benchmarking results for model V2. Top row: Temperature profile (left) and number densities of H and H2 (right). Middle row: Number densities 
of C+, C and CO (left) and surface brightnesses for [O I] at 63 jim and 146 jim; [C II] at 158 jim; and [C I] at 370 jim and 610 /an (right). Bottom row: 
local emissivities for [O I] 63 fim (left) and [C I] 370 fim (right). For this model we additionally compare 3D-PDR with the UCL_PDR11 code. 
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Figure 5. Benchmarking results for models V3 (left column) and V4 (right column). Top row shows temperature profiles, middle row shows number densities 
of H and H2 and bottom row shows number densities of C + , C and CO. 
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acting with a plane-parallel external radiation field (§43). In all 
these applications we use £ = levels of HEALPix refinement and 
#crit = 0.8 ~ 7r/4rad. The turbulent velocity is set to uturb = 
1 kms -1 . The dust temperature is fixed and set to Td ust = 20 K. 
The cosmic ray ionization rate is set to £ = 5 x 10~ 17 s _1 . In addi- 
tion we use Ichem = 8 iterations over chemistry at the beginning 
of each simulation and Ichem = 3 during each new iteration over 
thermal balance (see i]2.7t . 

We note that the applications presented here are simplified ex- 
amples which demonstrate however the capabilities of 3D-PDR in 
modeling any kind of density structure under the interaction of a 
UV radiation field. 



4.1 Interaction of a uniform-density spherical cloud with a 
plane-parallel radiation field 

In this application we consider a uniform-density spherically sym- 
metric cloud with a H-nucleus number density of n H = 10 3 cm~ 3 
and a radius of R — 5.15 pc. The cloud has therefore radial visual 
extinction of Av ~ 10 mag at its centre, assuming that the sur- 
face of the cloud is at Av = mag. A plane-parallel uniform UV 
radiation field of strength x — 10 Draines is impinging from one 
side. The density and UV field strength used for this application 
correspond to the parameters used in model VI of R07. 

We construct the sphere in the following way. Using HEALPix 
we create Mi = 3072 (level £ = 4) pixels uniformly distributed 
on the surface of the sphere and we take into account only those 
which lie on the hemisphere on which the UV field is imping- 
ing. Each one of these pixels defines the start of a line segment 
which penetrates the sphere; is parallel to the direction of the UV 
field; and consists of elements logarithmically distributed filling up 
the entire sphere. We use Na v = 60 elements per Av dex with 
Avmin = 10 mag which, as discussed in §3.21 ensures high res- 
olution along the direction of the UV field. Thus the total number 
of elements is approximately A4icm — 5.38 x 10 5 . 

Since this application is equivalent to Model VI in R07, in 
Fig(6]we compare our results with those of the benchmarked codes; 
we find in general a very good agreement, particularly with the 
UCLJPDR code; we note however that our temperature values at 
high Av are noisier primarily due to additional three-dimensional 
effects and due to our different iteration criterion which leads to 
less smoothing; we also find that our H/H2 transition occurs at 
slightly earlier Av- 

Figure [7] shows how the temperature varies when considering 
the limb and the equator of the sphere separately; as expected the 
temperature is slightly lower in the regions located around the limb 
of the sphere (as seen from the UV field) in comparison with the 
temperatures obtained in regions around the equator of the sphere. 
This is because the radiation field is impinging more radially in the 
latter than in the former leading to small differences in the attenu- 
ated field strengths. 



4.2 Interaction of a uniform-density spherical cloud with a 
two-component radiation field 

In this application we consider a uniform-density spherically sym- 
metric cloud with a H-nucleus number density of n H = 2 x 



10 cm and a radius of R = 2.58 pc. As in §4T| the cloud has a 
radial visual extinction of Av ~ 10 mag at its centre assuming that 
the entire surface is defined as Av = mag. A two-component 
UV radiation field was adopted. The first component corresponds 
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Figure 6. Results of Application 1 , in which we directly compare the one- 
dimensional codes with a fully three-dimensional calculation of 3D-PDR. 
From top to bottom: temperature profiles; number densities of H and H2; 
number densities of C + , C and CO. 
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Figure 7. Temperature profile of Application 1 obtained by 3D-PDR. Due 
to three-dimensional effects we find that the temperature of the limb is lower 
than that of the central region since in the latter case the UV field is imping- 
ing more radially. 



to a radial sampling field of strength Xiso = 120 Draines, and the 
second component corresponds to a plane-parallel radiation field of 
strength Xuni = 2 x 10 3 Draines impinging from one side. 

We constructed the sphere using a combination of two dif- 
ferent arrangements of elements. In the first arrangement we used 
HEALPix to create A/4 = 3072 (level 1 = 4) pixels uniformly 
distributed on the surface of the sphere. Each one of these pixels, 
along with the centre of the sphere, define ray segments which we 
constructed using Ma v = 20 elements logarithmically distributed 
per Ay dex and with Ay. min = 10~ 5 mag. Thus we created a 
sphere with approximately A4iem,i — 3.7 x 10 5 logarithmically 
radially spaced elements. In the second arrangement we create a 
sphere of A4leni,2 = 3 x 10 4 uniformly distributed elements. 
Therefore the resultant combined sphere consists of approximately 
A/" e iem = A/dem, i + A/dem, 2 — 4 x 10 5 elements and possesses 
both an approximately uniform distribution in its inner part and a 
logarithmic distribution in its outer parts, ensuring that the resolu- 
tion requirements described in i j3.2l are met. 

Figure [8] shows six different plots of our results. The plane- 
parallel radiation field is impinging from left to right. In the top 
and middle rows we plot the emission maps for [C II] at 158 /im 
(top left), [C I] 610 fim (top right), [O I] 63 (middle left) and 
CO (1-0) (middle right). At the bottom left we plot a cross section 
of the cloud, showing the gas temperature. We see that the temper- 
ature at the surface of the right-hand hemisphere is ~ 270 K due to 
the radial sampling component of the radiation field, and the tem- 
perature at the left-hand hemisphere reaches > 400 K due to the 
additional interaction of the plane-parallel component of the radi- 
ation field. The local undulations observed here do not correspond 
to local differences in cooling and heating but instead are a result 
of numerical noise introduced by the discretization in angle of the 
£ = choice of HEALPix rays. Although a selection of t > 
values would smear out these undulations, it would increase the 
computational cost without offering a significant improvement in 
the analysis. 

Overall, we see that the PDR is 'squeezed' at the left-hand 



side due to the plane-parallel radiation field. This is seen even for 
the CO (1-0) line which is embedded in the inner part. Converting 
to units of ergs -1 cm -2 sr -1 , we find that the strongest coolant 
is the [O I] 63 /im line, which on average is ~ 1.7 times stronger 
than the second coolant, [C II] 158 jLtm. At the bottom right we 
show an RGB composite image of the CO (1-0) (red), [C I] (green) 
and [C II] (blue) emission maps. Here, we observe a well-defined 
stratification of species. 

In this application we additionally explore how the calcula- 
tions of the one-dimensional treatment of PDRs diverge from the 
corresponding calculations when a fully three-dimensional treat- 
ment is taken into consideration. To do this we perform a one- 
dimensional calculation of a PDR which has the same parameters 
of those explored above (i.e. n = 2 x 10 3 cm -3 , Na v = 20 el- 
ements logarithmically distributed per Ay dex and with Ay iXOin = 
10~ 5 mag and Ay. n ,„ = 10 mag, and \ = 2120 Draines field 
strength impinging from one side). We then compare the attenua- 
tion of the UV field strengths and the number densities of C + , C, 
and CO of this run and the corresponding values taken from two 
different radial directions from the spherical cloud. These direc- 
tions are shown in the RGB composite image of Fig[8] The dotted 
line (Equator) is parallel to the direction of the plane-parallel UV 
field and the dashed line (Diagonal) is not. In particular the Equator 
corresponds to the HEALPix ray ID=6 of I = and the Diagonal 
to the ray ID=9, supposing that a 12-ray structure is emanated from 
the centre of the sphere. While the radial visual extinction of the 
sphere is quite high, we neglect the contribution of the radial sam- 
pling of the UV radiation impinging from the opposite side of the 
cloud. 

Figure [9] shows the attenuation of the UV field (upper panel) 
and the number densities of C + ,C and CO (lower panel) for the 
one-dimensional calculation (solid line), along the Equatorial di- 
rection (dotted line), and along the Diagonal direction (dashed 
line). Although the unattenuated field strength at the surface of 
the PDR is the same in all cases, the difference due to the three- 
dimensional structure has an impact in the attenuation of the UV 
field as seen from the Diagonal direction. On the other hand, due to 
the symmetry obtained, the Equatorial direction is in a very good 
agreement with the one-dimensional calculation. We therefore find 
a difference in the distribution of the C + , C and CO abundances 
depending on the direction along which we perform an observation 
even in the present simplified case. 



4.3 Interaction of a cometary globule with a plane-parallel 
radiation field 

In this application we explore the capability of the code to simulate 
PDRs with arbitrary density distributions. To do this we used as ini- 
tial conditions a snapshot taken from a Smoothed Pa r ticle H ydro- 
dynamic simulation presented in §4.4 of Bisbas et al. U2009I) That 
SPH simulation (using the SEREN code; Hubb er et alj|201lh ex- 
amined the interaction between an initially uniform density clump 
with UV radiation emitted spherically from an exciting sou rce (by 
invoking the on-the-spot approximation; IOsterbrock|[l974h which 
wa s placed outside an d far away from the clump (for full details 
sec iBisbas et alj|2009b. This interaction, referred to as " radiation 
driven implosion" jS andford, Whita ker. & Kleinl 1 19821 ; iBertoldil 
1989; Lefloch & Lazarefn ll994). drives a strong shock front into 
the inner part of the clump, creating a morphological structure rem- 
iniscent of the cometary globules obser ved in the interior of io nized 
regions (such as in the Helix Ne bula; |Matsuura et alJl200S ), and 
which may trigger star formation ( Kessel-Devnet & Burkert 2003; 
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Figure 8. Results of Application 2, in which we simulate the interaction of a spherically symmetric cloud as it interacts with a UV field consisting of a 
radial sampling and a plane-parallel component. The top four plots show emission maps for [C II] 158 fim (top left), [C I] 610 (im (top right), [O I] 63 fj,m 
(middle left) and CO (1-0) (middle right). The colour bars are in units of (K km s _1 ). The bottom left image shows a cross section of the gas temperature 
(K). The bottom right image shows an RGB composite image of the emission maps of CO (1-0) (red), [C I] (green) and [C II] (blue). The values on the colour 
bar correspond to the [C II] emission. RGB colour bar ratios of 5:1:10 for CO(1-0):[C I]:[C II]. The white dotted line in the RGB image corresponds to the 
direction along the "Equator" whereas the dashed is along the "Diagonal" direction (HEALPix ray ID=6,9 of the NESTED numbering scheme - see i|4.2l for 
the relevant discussion). 
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Figure 9. Divergence of calculations between one-dimensional and fully 
three-dimensional calculations of the PDR described in i|4.2l The top panel 
shows the attenuation of the UV field strength (in units of the Draine field) 
and the bottom the corresponding number densities of C + , C, and CO (in 
cm -3 ). The solid line (ID calc.) is the one-dimensional calculation and 
the dotted (Equator) & dashed (Diagonal) lines are the directions shown 
in the RGB composite image on bottom right of Fig[8] When effects of 
the three-dimensional structure of the PDR are taken into consideration, 
discrepancies between the calculations appear. 



Gritschneder et alj 120091 ; iBisbas et al] 1201 ll ; lHaworth & Harries! 
2012b . 

We take a snapshot at t — 0.12 Myr. A cross-section plot (at 
z = 0) of the density structure of the clump at t — 0.12 Myr is 
shown at the bottom left of Fig[l0] where the ionizing radiation is 
impinging from bottom to top. At this time the cloud has attained a 
V-shape structure which contains an approximately ellipsoidal core 
of density n H ~ 10 5 — 10 6 cm" 3 located at its tip and directly ex- 
posed to the ionizing radiation. Although in the SPH simulation the 
radiation is emitted spherically symmetrically from the distant ex- 
citing source, the angular size of the clump is quite small (~ 6°) 
and we therefore consider a plane-parallel radiation field here. The 



UV photon flux impinging is $ ~ 2.18 x 10 9 cm -2 s _1 corre- 
sponding to a field strength of approximately \ = 30 Draines. 

For the purpose of this application and for the sake of com- 
putational speed, we only perform calculations up to Av )BB = 
4 mag of visual extinction. For Av > 4 mag and for the radiation 
field strength and density considered for this application, the cloud 
cannot be treated as a PDR anymore as it has reached dark cloud 
conditions, i.e the radiation field does not penetrate any longer and 
hence does not have any effect on the chemistry. Instead, the lat- 
ter will be mainly dominated by cosmic ray induced reactions (in- 
dependent of optical depth). For a proper treatment of dark cloud 
chemistry, depletion on to dust grains should be taken into consider- 
ation. The temperature has reached equilibrium values of ~ 10 K; 
CO lines are mainly saturated and hence not contributing much to 
the cooling, and no strong source of heating contributes to increas- 
ing the temperature. To avoid calculations in this dark molecular 
region, we use the following technique. For a random element, p, if 
the magnitude of Av of the HEALPix ray with the highest attenu- 
ated flux exceeds a user-defined threshold (here Ai/ !m „ = 4), then 
p is not considered as a PDR and calculations are omitted. Here, 
we also omit all regions with n H ^ 100 cm -3 as these belong to 
the inner part of the H II region. 

Figure [10] shows the results of our calculations. The top four 
frames show the emission maps for [C II] 158 fjxa (top left), [C I] 
610 Atm (top right), [O I] 63 pm (middle left) and CO (1-0) (middle 
right). In these maps we plot only the PDR (the contribution due to 
the dark molecular component is excluded). At the bottom right we 
show an RGB composite image of the emission maps of CO (1-0) 
(red), [C I] (green) and [C II] (blue). 

The density structure shows some symmetry with no abrupt 
gradients, i.e. the density changes smoothly and there are no sharp 
density enhancements. From the emission maps we see that the 
species in the PDR are distributed smoothly and follow the density 
profile. The weakest emission is produced by [O I] 63 fim and the 
strongest by CO (1-0) implying that the molecular gas dominates 
over the atomic contribution. However, considering the transition 
frequencies for these maps, we find that the [O I] 63 fim and the 
[C II] 158 fim lines are the dominant coolants with the first being 
somewhat (< 10%) stronger. In addition, the thickness of the PDR 
is very small comparing with the dimensions of the cometary glob- 
ule. This is observed in the RGB composite image. Although there 
is a stratification of the species with [C II] emitted from the outer- 
most part and CO (1-0) from the innermost part of the globule as 
expected, the transition between the species occurs in quite a thin 
layer. Since the colouring of these species overlap, we observe a 
bright white rim around the whole cometary globule. 



5 CONCLUSIONS 

We have presented 3D-PDR, a numerical code for simulating three- 
dimensional PDRs of arbitrary density distribution. The code uses 
a ray-tracing scheme based on HEALPix in order to calculate the 
column densities and the escape probability in every direction for 
every element within the cloud. We adopted a reduced chemical 
network of 33 species and 320 chemical reactions. Through an it- 
erative process the code calculates the cooling and heating rates for 
every cloud element adjusting the gas temperature at each iteration 
in an attempt to balance the heating and cooling rates. The code 
terminates once the difference between the heating and cooling is 
negligible i.e. when the PDR has obtained thermal balance. 

We tested the ray-tracing scheme in calculating the column 
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Figure 10. Results of Application 3, in which we use as initial conditions a snapshot from a Smoothed Particle Hydrodynamics simulation of a cloud as it 
interacts with UV radiation, undergoing radiation driven implosion. The radiation is impinging from bottom to top. The top four plots show emission maps 
for [C II] 158 /im (top left), [C I] 610 fim (top right), [O I] 63 /im (middle left) and CO (1-0) (middle right). The colour bars are in units of (K km s _1 ). 
The bottom left image shows a cross section of the number density profile of the cloud. The bottom right image shows an RGB composite of emission 
maps of CO (1-0) (red), [C I] (green) and [C II] (blue). The values on the colour bar correspond to the [C II] emission. RGB colour bar ratios of 8:1:2 for 
CO(1-0):[C l]:[C II]. 
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density of a particular element against its known analytical expres- 
sion and we found very good agreement. We have also explored the 
spatial resolution requirements for simulating a PDR and we found 
that our code resolves one-dimensional uniform-density PDR if it is 
constructed using Ma v = 20 elements logarithmically distributed 
per Av dex. 

Furthermore, we repeated the benchmarking tests presented in 
R07 and we compared our results with the ones obtained by other 
one-dimensional PDR codes. Overall we find very good agree- 
ment between the one-dimensional codes and 3D-PDR. In addi- 
tion, we explored the capabilities of our code in simulating three- 
dimensional structures exposed in one- or two- component UV ra- 
diation fields. In particular: 

• We examined the interaction of a uniform-density spherical 
cloud with a plane-parallel radiation field in which the values of 
density and field strength were identical to those of model VI in 
R07. We found very good agreement between 3D-PDR with the 
one-dimensional codes and in addition we observed at low Ay 
cooler temperatures in the limb of the sphere in contrast with the 
higher temperatures in the equatorial regions; an effect directly re- 
lated to the three-dimensional treatment in our case. 

• We examined the interaction of a uniform-density spherical 
cloud with a two-component radiation field, consisting of an ra- 
dial sampling field and a plane-parallel field. We explored the dif- 
ferences in results obtained when a fully three-dimensional treat- 
ment of the PDR is taken into consideration in contrast with a one- 
dimensional simplification. We found that the results differ accord- 
ing to the direction at which observations are performed. 

• We examined the interaction of a cometary globule with a 
plane-parallel radiation field, where we considered as initial condi- 
tions a snapshot taken directly from an SPH simulation. We found 
that the PDR location follows the density profile of the globule, i.e. 
the abundances of species change in agreement with the density 
structure. We also found that the thickness of PDR is quite small 
in comparison with the overall size of the globule using composite 
RGB emission maps. This application showed also the capability 
of our code to model any type of density and particle distribution. 

The coupling of 3D-PDR with MOCASSIN will be presented 
in a forthcoming paper. The integrated code should make feasible 
a realistic treatment of three-dimensional H II/PDR complexes in- 
cluding a detailed treatment of SEDs, thus offering a powerful tool 
in studying such structures with arbitrary density distributions and 
multiple exciting sources. 
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The code starts by reading the inputs for the model, includ- 
ing the density structure of the cloud, the initial abundances, and 
the physical parameters describing the environment; it builds the 
evaluation points using a HEALPix based ray-tracing scheme (see 
§2.21 ; it calculates the total column density to the cloud surface for 
each spatial element within the cloud; and computes the attenua- 
tion of the user-defined UV radiation field (see E|2.3t . Currently, the 
user is able to choose between three simplified types of UV field: 
UN I -directional (plane-parallel), isotropic, and PoiNT source, or 
any combination thereof. 

An initial "guess" temperature is assigned to each cloud ele- 
ment based on the UV field strength at that point (see j|2.8) . The 
chemical reaction rates and resulting abundances and column den- 
sities for each species are then calculated, and this cycle is repeated 
Ichem times to reach convergence ( §2.7\ . LTE level populations are 
then determined for each coolant species as initial guesses for the 
radiative transfer calculation. 3D-PDR applies a three-dimensional 
escape probability method (see £]2.4| > to treat the line transfer, up- 
dating the level populations accordingly and iterating to obtain 
level population convergence (judged according to a user-defined 
tolerance parameter, as described in S |2.8t . 

Once converged, the total cooling and heating rates are com- 
puted from the sum of the individual contributions (see i]2.5l and 
£]2.6I > and compared to determine if thermal balance has been 
reached ( i]2.8l >. If not, new temperatures are assigned to each cloud 
element and the iterative search for thermal balance continues. 
Once thermal balance has been reached at all cloud elements, the 
code writes the outputs and terminates. 



APPENDIX A: FLOWCHART OF 3D-PDR 

Figure IaTI shows a flowchart of the computational scheme used in 
3D-PDR. Each solid box corresponds to a DO-loop over all ele- 
ments within the cloud and the dashed box corresponds to iteration 
on the chemistry. 
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Figure Al. Flowchart of 3D-PDR 



